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ABSTRACT 

Motivation: Gene transcription requires the orchestrated binding 
of various proteins to the promoter of a gene. The binding times 
and binding order of proteins allow to draw conclusions about the 
proteins' exact function in the recruitment process. Time-resolved 
Chip experiments are being used to analyze the order of protein 
binding for these processes. However, these ChIP signals do not 
represent the exact protein binding patterns. 

Results: We show that for promoter complexes that follow sequential 
recruitment dynamics the ChIP signal can be understood as a 
convoluted signal and propose the application of deconvolution 
methods to recover the protein binding patterns from experimental 
Chip time-series. We analyze the suitability of four deconvolution 
methods: two non-blind deconvolution methods, Wiener deconvolution 
and Lucy-Richardson deconvolution, and two blind deconvolution 
methods, blind Lucy-Richardson deconvolution and binary blind 
deconvolution. We apply these methods to infer the protein binding 
pattern from ChIP time-series for the pS2 gene. 
Contact: rudolf.hanel@meduniwien.ac.atj 



An essential step in gene expression is the initiation of transcription. 
This process requires the orchestrated recruitment of regulatory 
proteins to the gene promoter. This leads to the formation of the 
transcriptional machinery that transcribes a gene from DNA to 
RNA. Transcription factors assemble on the promoter site, forming 
sequences of protein complexes on the promoter Eventually, the 
protein complex attracts the RNA polymerase that transcribes 
the gene and finally the promoter is cleared again. During this 
process the modification of epigenetic marks on the DNA and 
histones were shown t o be essential for the initiation of transcription 
( JMetivieretailiaOOel) . 

For various promoters the proteins participating in the complex 
formation have been identified, yet, determining the exact order and 
timing of the recruitment events is still a non-trivial task. 

Extensive ChIP experiments have been used to examine the 
binding order of prote i ns and the modifi cation of epigentic marks 
JMetivier et~an. l2006l : JLee et al.L |2005|) . Since a large number 
of cells (~ 10^) a re required to perform ChIP experiments 
( JMetivier et aliuOOq) . such experiments are performed on initially 
synchronized cell populations with a cleared promoter site. 
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Fig. 1. Visualization of the recuritment network for sequential recrutiment. 
The recruitment matrix aji circularly links the recmtiment states and each 
state has only one successor state. 



Usually Chip experiments are analyzed heuristically by applying 
prior knowledge of protein interactions to interpret the form of the 
Chip signal. 

In ? the sequence of maxima and minima in ChIP time-series was 
used to predict the structure of the dominant negative feedback loop 
that governs the recruit ment dynamics. 

In lHanel et alj (2012) ChIP data was analyzed by representing the 
recruitment process as a regulatory network. The proteins involved 
in the recruitment process are represented by nodes in this network. 
By linearizing the regulatory dynamics, one can determine how the 
binding of a protein affects the affinity of other proteins to bind. 

By assuming a circular recruitment process that is traversed 
stochastically in each cell, ? reproduced the binding pattern by 
a least-square fit of a ChIP signal against simulation data of the 
recruitment process. This method results in binding patterns that 
show multiple binding times for most proteins for the pS2 promoter. 
These binding patterns have been interpreted as stochastic binding 
events that are seen for example with other experimental methods 
like GFP (^reen fluorescent proteins). 

In Scholling et alj <2013h we generalized this recruitment model 
to address the question of whether the binding order in a recruitment 
process is deterministic (sequential recruitment) or stochastic 
(probabilistic recruitment). In this model the recruitment process 
is represented by a walk on a network of recuritment states. In the 
case of sequential recruitment each recruitment state has only one 
possible successor state (compare Fig. [T) whereas in the case of 
probabilistic recruitment the topology of the network can look much 
more complex. We showed that in a cell population the occupation 



) Oxford University Press 2013. 
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Pj (t) of a recruitment state j at time t is given by 



^Pj(*) = 5Z ("j* ~ ak5j,k)Pk{t) 



(1) 



where a^fc characterizes the transition rate from the recuritment 
states k to state j, 5j_fe is the Kronecker symbol and a^ = X^i 'isfc- 
If a,fc = 0, a t r ansitio n from state k to j is impossible. Following 
IScholling et alj ( 120131) the resulting ChIP signal Cm (t) for a protein 
m in a circular sequential recruitment process with frequency ^ can 
be predicted by 



.(t) 



Sm(x - flt)Po{x)dx 



(2) 



in the limit of a large number of recruitment states (TV — >■ cxd), where 
the periodic function Sm{x) = Sm{x + 1) denotes the protein 
binding pattern of a protein (Sm{x) = 1 if protein m is bound in 
state x and otherwise), x = k/N and Po{x) = PxN(t = 0). The 
assumption of large N can be justified in the context of the pS2 gene 
IScholling et al., 2013; Metivier et al., 2003). 

Consequently, the ChIP signal Cm (t) does not directly represent 
the binding pattern S but it is "blurred" by Pq. 

Based on this observation, we analyze methods to recover the 
protein binding pattern from ChIP time-series by deconvolution 
in this paper. Deconvolution methods attempt to recover the 
original, un-bluiTed source signal from a measured blurred 
signal. The problem of deconvolving signals is fou nd in various 
fields of research, r a nging f rom imag e proc essing (IFergus et all 
l2006l: ISchulzetali 1 19971 : ICamionf . Il976h over engineer ing 
I McDonald et al. . 2012 ) a nd spectroscopy JNadler et al.iri98 9l) to 



biologv toown et al.Ll200ilLun et al.Ll2009h . In the context of ChIP 
time-series, we assume that the synchronization of a cell population 
is imperfect, e.g. due to variation in th e time required by each cell 
to translocate proteins into the nucleus ( lAshall et al.ll2009l) . 

Deconvolution can be performed using various mathematical 
approaches, which are based on different assumptions on the 
source signal and/or the blurring process. In this paper we 
discuss the application of four deconvolution methods to recover 
the protein binding pattern from ChIP time-series: Wiener 
deconvolution, Lucy-Richardson deconvolution, blind Lucy- 
Richardon deconvolution and binary blind deconvolution. 

We demonstrate the application of these four deconvolution 
metho ds on ChIP time-series data for the pS2 gene ( Met ivier et al.i, 
120030 and discuss the results of the different methods. 



1= 50 


____..^ V 


_..^^^"V 


^--^1 


5D 


100 




S'i 5° 










50 


100 


15 


S| 50 


-~. -.--^ "■^ 


__^ ^ V. 


____/' 




50 


100 


15( 


1" " 


y. V 


_y^^'-^ 


/- " 




50 


100 




Si" 




~"""""'""\. 






50 


100 


15 


£| 50 


y v^ 


^^--x 


_^^ 




50 


100 


15( 


^ 








50 


100 




is » 










50 


100 


15 


e| 50 






""\ /■ 


» 


50 


100 


15( 



Fig. 2. Chip time-series for the pS2 promoter: The time-series show the 
Chip signals (normalized) for five proteins and modifications (methylation 
and acetylation) of two histones. Each time-series covers the first, 
unproductive cycle and three productive cycles (data extracted from 
lMaiviereta"i]J20O3h : ?). 



state (compare Fig. [T}. We showed that the probability of a cell to 
be in state j is then given by Eq. (|T). Since the number of cells 
in a Chip experiment is large, pk{t) does not only represent the 
probability for a single cell to be in state k but also estimates the 
probability of occupation of a state within the whole cell population 
at a given time. 

The infoiTnation of whether a protein m is present in a state k is 
represented by a matrix Smk '■ If the protein is bound in state k then 
Smk = 1 and Smk = otherwise. By calculating the probability 
Pk (t) of a cell to be in state k at time i, the resulting ChIP time-series 
Cm (t) can be predicted by 



Cm{t) = 22 ^"^kPkit) . 



(4) 



1 METHOD 

In lSchoUing et al.l J2013r) we represented the recruitment process as 
a walk on a network of recruitment states. Each cell traverses this 
network, one state after another. The waiting times r between the 
state transitions from i to j are considered as a probabilistic event 
drawn from an exponential distribution 



Pji{r) = ajiexp(-ajir), 



(3) 



where the aji deteiTnine the transition rates. In particular if no 
transition i -^ j is possible, then aji — 0. In the case of 
sequential recruitment each state has only one possible successor 



For a large number of recruitment states (N -^ oo) one obtains that 
the Chip signal for a sequential recruitment process is a convolution 
of the protein binding pattern Sm{x) with the initial distribution of 
states Po{x), where is constant over time (Po{x) = pxN{t — 0), 
see Eq. l|2j)- Thus the ChIP signal does not directly represent the 
binding pattern but is "bluiTed" by Po{x). 

Po{x) can be identified as the convolution kernel, the Sm{x) as 
the deconvoluted signal and the Cm(i)'s as the convoluted signals. 
One can therefore apply deconvolution methods to recover the 
protein binding pattern from the ChIP signal. 

Deconvolution methods are frequently used in image processing 
where the task is to reconstruct the original image from a blurred 
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Fig. 3. Results for Wiener deconvolution (5% signal-to-noise ratio, V row: deconvolution kernel). P' column: ChIP time-series for five proteins and 
modifications (methylation and acetylation) of histones H3 and H4. 2"'' - 6* column: protein binding pattern (nomialized) for the deconvolution in Wiener 
space with Gaussian kernel with standard deviations cr = {1 min, 2 min, 3 min, 5 min, 7.5 min}. 



image. From a mathematical point of view, the task is to recover the 
unknown, un-blurred function g{x) = Sm (i) from a known, blurred 
function f{x) = Cm{t) where h{x) describes the blurring process. 
The convolution kernel h{x) is unknown in many applications. To 
overcome this dilemma, one either has to assume a certain kernel 
function or use methods that are able to estimate the kernel h{x) 
from the blurred image with the help of certain assum ptions. The 
latter methods are calle d blind deconvolution methods JFish et aI.L 
ll995l:lLevinetal.Ll2009l) . 

In this paper we assume that the kernel is given by a Gaussian 
distribution, meaning that some cells are already ahead in the 
re cruitment process wh ereas others lag behind. Bas ed on the results 
of lMetivier et al.l l l2003h and lSchoUing et alj ( 1201 3h . we assume that 
the recruitment process is sequential and the kernel is constant in 
time, i.e. the de-synchronization happened prior to the considered 
time interval that is used for deconvolution. The latter can be 
justified by the fact that no obvious decay of signals is de tectable 
in the ChIP time-series we consider (iScholling et al.Ll2013h . Thus, 
no further (de-)synchronization of the cell population occurs during 
this time interval. 

Note that the analysis of recruitment processes that are non- 
sequential, i.e. these processes are not traversed in a cyclic order but 
contain "shortcuts" or alternative recruitment pathways, is beyond 
the scope of this paper. 

An important property of ChIP assays is their sensitivity: the 
unit of Cm{t) is expressed in "percentage of input", meaning the 
relative abundance of a specific DNA sequence using a protein- 
specific antibody compared to a non-specific antibody control. This 
approach leaves one degree of freedom in amplitude of the ChIP 
time-series since different antibodies can have different affinities 
and thus cause ChIP signals that differ in amplitude. Consequently, 



this property also transfers to the amplitude of the protein binding 
pattern Sm (t)- For this reason, we normalize the amplitude of ChIP 
time-series Cm{t) and the protein binding pattern Sm{t)- 
In this paper we consider four deconvolution methods: 

• Wiener deconvolution (?) is a deconvolution method that works 
in the frequency domain. In contrast to naive deconvolution 
in Fourier space, it suppresses amplification of errors at 
frequencies with a low signal-to-noise ratio. 

• The Lucy-Richardson (LR) deconvolution is a method based 
on Bayes' theorem. It maximizes the likelihood of the restored 
signal to b e the tr ue one by using the expectation-maximization 
algorithm iILucvI 1974). 

• Blind Lucy-Richardson deconvolution tries to recover the 
convolution kernel by exploiting the commutative property 
of the convolution operation: With an initial guess of the 
kernel, the LR scheme is applied iteratively first to calculate 
a preliminary deconvolution of the convoluted signal and then 
to calc ulate a new kernel from th is preliminary deconvolution 
result CBig gs and Andrewslll997h . 

• Binary blind deconvolution JLaml 120071) applies a binary 
constraint on the deconvoluted signal (g{x) G {0, 1}). The 
deconvolution can then be cast into a convex optimization 
problem, which can be solved by standard methods for convex 
optimization. 

In contrast to the first two deconvolution methods the latter 
methods are blind deconvolution methods: These algorithms try to 
predict a convolution kernel h(x) instead of using a given kernel as 
input. The problem of blind deconvolution is ill-posed in general 
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Fig. 4. Results for Lucy-Richardson deconvolution (same setup as in Fig. (3))- For Gaussian kernels with small standard deviation the binding patterns do not 
differ strongly from the corresponding ChIP time-series. The protein binding patterns for larger Gaussian kernels show distinct peaks for the binding events. 
The protein binding patterns for TBP and Met-H4 show plateaus. 



JLevin et all 120091) : Many pairs {h{x),f{x)} of kernel and un- 
blurred function exists that can reproduce the observed, bluiTed 
signal f{x). Thus additional assumptions on the kernel and/or on 
the un-blurred signal are necessary. A detailed description of the 
deconvolution methods can be found in the Supporting Information. 
We investigate the suitability of these four deconvolution methods 
to recover the protein binding pattern from ChIP signals for the pS2 
gene in the next section. 



2 RESULTS 

We apply the four deconvolution methods to ChIP time-series for 
the promoter of the p S2 gene in human MCF-7 breast cancer cells 
<M6tivier et aI.L 120030 . Transcription of the pS2 gene is induced by 
the human estrogen receptor hERa. Here we consider the binding 
of five proteins (hERa, HDAC, PolII, TBP, TRIPl) and the histone 
modifications (acetylation and methylation) of two histones {H3 and 
H4). The time-series of the ChIP data in Fig. l|2j covers a time 
interval of 145 minutes with an experimental temporal resolution 
of approximately 5 minutes and exhibits oscillatory dynamics. We 
resampled the data to a time axis with Ai = 15sec to account for 
non-equidistributed time points. Since the first transcription cycle 
is unproductive, we only consider the time span 35~145min of the 
experiment. A striking featu r e of t he ChIP signals in the dataset 
extracted from lMetivier et al.lcOOSh is that two of the signals {Met- 
H4 and TBP) show distinct dynamics from the other five signals: 
Instead of three oscillations, only two oscillations can be observed. 
Note that exclusion of these two signals from our analysis does not 
alter the deconvolution results of the other signals. 



2.1 Wiener Deconvolution 

The results for Wiener deconvolution are shown in Fig. (|3] for 
Gaussian kernels with standard deviations g G {1, 2, 3, 5, 7.5} min. 
We estimate the signal-to-noise ratio to be 5%. A signal-to-noise 
ratio of 10% only alters the deconvolution results slightly (not 
shown here). For kernels with low standard deviation (cr < 
2 min) the deconvolution results barely differ from the ChIP time- 
series. For larger standard deviations the deconvolution results in 
protein binding patterns exhibit negative amplitudes. This violates 
the constraint that concentrations can only be non-negative and 
therefore have to be discarded. The results of Wiener deconvolution 
show that this deconvolution method only results in valid binding 
patterns for Gaussian kernels with small standard deviation (a < 
2 min) and these binding patterns show only marginal differences 
compared to the original ChIP signals. 



2.2 Lucy-Richardson Deconvolution 

The results of the Lucy-Richardson deconvolution method are 
shown in Fig. Q. The deconvolved signals for Gaussian kernels 
with a small standard deviation {a < 10 min) do not differ strongly 
from the ChIP signals. Artifacts in the form of large peaks are 
visible at the end of the protein binding patterns at f = 145 min 
for a G {1 min, 2 min}. However, for a Gaussian kernel with 
large standard deviation (a = 10 min) the deconvolution results 
in distinct peaks for the protein binding patterns. Inbetween these 
two peaks the protein binding patterns decrease to zero (with the 
exceptions of TBP and Met-H4 as we discuss below). 
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Fig. 5. Results for Lucy-Richardson deconvolution with the Gaussian kernel a = 9.69min) ( 1"' row: deconvolution kernel). P' column: ChIP time-series 
for five proteins and modifications (methylation and acetylation) of histones H3 and H4. 2"''-6"' column: protein binding pattern (normalized) for the Lucy- 
Richardson deconvolution with Gaussian kernel with standard deviations a = 9.69 min. 3"* row: result of re-convolving the deconvolved signal with the 
Gaussian kernel (compare to P' column). 



Since convolution blurs the signal, we searched for the Gaussian 
kernel that results in the sparsest deconvolved signal. This Gaussian 
kernel was found at cr = 9.69 min and is shown in Figs. l5]|6j. 

The signals that show approximately three periods in the 
Chip time-series, also show three distinct binding times in the 
deconvolution result. The ChIP signals for methylation of H4 and 
binding of TBP show conceptually different dynamics than the 
Chip signals of the other proteins (two oscillations instead of 
three oscillations). As a consequence, their binding patterns also 
show different characteristics than the binding patterns of the other 
signals: The first and the second binding events are not very distinct 
in time. Inbetween these two peaks the binding pattern exhibits 
plateaus in the affinity of about 0.2 and 0.15 (arbitrary units) for 
TBP and Met-H4, respectively. 



2.3 Blind Lucy-Richardson Deconvolution 

The results for the blind Lucy-Richardson deconvolution are 
shown in Fig. Q. Note that blind Lucy-Richardson deconvolution 
starts from an initial guess of the kernel and estimates the un-blurred 
signal and the convolution kernel with the maximal likelihood. 
We started the deconvolution with Gaussian kernels. The resulting 
kernels are shown in the first row of Fig. (|7) indicated in red. The 
deconvolved signals are very similar to the original ChIP signal 
for all kernels. Interestingly, for initial Gaussian kernels with large 
standard deviation (a > 15 min) the initially broad kernels are 
always driven towards sharply peaked distributions, which explains 
the small effect yielded by this method. 



2.4 Binary Blind Deconvolution 



The deconvolution algorithm proposed bv lLaml(l2007h performed 
worst: The algorithm failed to converge for the ChIP time-series for 
initial Gaussian kernels. For this reason, we can not present any 
protein binding patterns for this deconvolution method. 



3 DISCUSSION 

Three of the four deconvolution algorithms we tested were able 
to deconvolve the protein binding patterns from the ChIP signals. 
However, Wiener deconvolution only producted valid binding 
patterns for sharp convolution kernels (small a) and also blind Lucy- 
Richardson deconvolution resulted in sharp convolution kernels. 

As one may expect deconvolution with these kernels does not 
show substantial differences to the original ChIP signals. This can 
be explained by the fact that these Gaussian kernels are very near 
to the Dirac delta function 5{x). Note that convolution of a signal 
with the Dirac delta function results in the very same signal again 
{{f®5)ix)^f{x)). 

Non-blind Lucy Richardson deconvolution resulted in protein 
binding patterns that did substantially differ from the ChIP signals. 
Regarding the results for the kernel that reproduces the sparsest 
admissible protein binding pattern (a — 9.69 min), the first events 
that occur in the recruitment process are binding of TBP and 
methylation of histone H4. However, these events are not indicated 
very sharply in the deconvolved signal and thus their binding times 
are less reliable. We defer the discussion of the deconvolution results 
for these events to the end of this section. 
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The protein binding patterns for the other four binding and histone 
modification events are indicated more distinctly than for Met-H4 
and TBP. According to the deconvolutions in Fig. (|6j, the event 
that occurs first is the binding of hER(W This result is consistent 
with prior biological knowledge, since the transcription of pS2 is 
known to be initiated by the binding of this estrogen receptor. Next 
in the recruitment process the histone H3 is methylized, which is 
indicated by sharp peaks in each of the three recruit ment cycles. 
CBP is known to be capable of methylizing this histone JWang et all 
l200lh . However, CBP is not included in the present experimental 
data and thus we cannot check whether CBP is accountable for this 
methylation event. 

Deacetylation events of H3 and H4 follow according to the results 
of LR deconvolution. These events end just w hen the binding pattern 
of RNA polymerase II peaks. Interestingly, iMetivier et al.l yOOo) 
defines a transcriptionally engaged pS2 promoter by the presence 
of Met-H4 and Ac-H3. In the results of LR deconvolution the 
methylation time of H4 is not indicated very precisely, but the 
time of deacetylation for H3 matches with the binding of RNA 
polymerase II. 

The re-convoluted signals of HDAC and TRIPl exhibit the largest 
errors in reconstructing the initial ChIP time-series (Fig. [5]l and 
regarding to our results they bind very late in the recruitment 
process. It is known that these proteins are involved in promoter 
clearance ( Metivier et al.L l2003n which is in agreement with the 
late binding time predicted by our results. However, the predicted 
binding time of HDAC (histone deacetylase) does not coincide with 
the deacetylation of H3 and H4. 

Thus we can propose two hypotheses: Either this means that 
HDAC does not require to bind to the protein complex to deacetylize 
the histones or HDAC is bound while histone deacetylation but the 
Lucy-Richardson prediction of the beginning of the HDAC binding 
time is incorrect. 

In this case the binding of HDAC better coincides with the 
deacetylation events for kernels with small standard deviation 
(cr < 2 min) e.g. as predicted by blind Lucy-Richardson. Thus the 
second possibility would indicate that cell sample used in the ChIP 
experiment is strongly synchronized. 

We return to the discussion of Met-H4 and TBP, which 
show distinct dynamics in the ChIP signals: instead of three 
oscillations like the other time-series they only show two 
oscillations. IMetivier et alj ( 120030 argued that detachment of 
TBP and demethylation of H4 do not occur in each productive 
recruitment cycle. Instead these events occur only in every second 
cycle. Besides the strict alternation of even and odd cycles, also 
probabili stic alternation is a pos sible interpretation of these binding 
patterns JScholling et al.L l20I3h . In this context, one can interpret 
the plateaus between the two peeks for TBP and Met-H4 in Figs. ([H 
|6j as signals only from those cells in which TBP does not detach 
and H4 is not de-methylized after the first productive cycle. 




Fig. 6. Results for Lucy-Richardson deconvolution with the Gaussian kernel 
a = 9.69 min (shows the same data as 2"'' column of Fig. |5): blue areas 
indicate low binding affinity, red areas indicate high binding affinities. 



4 CONCLUSIONS 

In the present paper, we have analyzed the applicability of 
different deconvolution methods to recover the protein binding 
pattern from ChIP time-series with an imperfectly synchronized 
cell population. Our analysis showed that Lucy-Richardson 
deconvolution can be used to improve inference on protein binding 
pattern from ChIP time-series. The presented results allow a 
more precise inteipretation of the ChIP data in forms of the 
function proteins play in the recruiting protein complex and 
the recruitment process in general. We considered two non- 
blind deconvolution methods, Wiener deconvolution and Lucy- 
Richardson deconvolution, as well as two blind deconvolution 
methods, blind Lucy-Richardson deconvolution and binary blind 
deconvolution. The applicability of these methods was tested on 
Chip data for the pS2 gene. Binary blind deconvolution failed 
to perform deconvolution. Wiener deconvolution and blind Lucy- 
Richardson deconvolution resulted in protein binding patterns 
that did not differ substantially from the original ChIP signals. 
For Wiener deconvolution broader convolution kernels violate the 
positivity constraint of the deconvoluted signal. However, the 
Lucy-Richardson deconvolution method has proofed its potential 
to recover protein binding patterns with distinct binding times. 
The binding pattern of HDAC indicates a late binding time which 
corresponds to its clearance function in the recruitment process but 
does not coincide with the deacetylation events of H3 and H4. This 
can either be interpreted as an inaccurate prediction of the proteins' 
binding pattern for this deconvolution method or it means that 
HDAC does not require to bind to the promoter to deacetylize the 
histones. The reconstruction of the binding patterns succeeded for 
all other proteins. Their timing is in agreement with prior biological 
knowledge and can help to investigate the exact function of the 
proteins in the recruitment process. 

Funding: This work has been supported by the Forum Integrativ- 
medizin an initiative of the Hilde Umdasch Privatstiftung. 



' For the iirst productive cycle one can question whether the binding of 
hERa or the methylation of histone H4 occurs first. But in the second and 



third cycle a clear distinction is possible. 
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